Quantum initial condition sampling for linearized density matrix dynamics: 
Vibrational pure dephasing of iodine in krypton matrices 
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This paper reviews the linearized path integral approach for computing time dependent properties 
of systems that can be approximated using a mixed quantum-classical description. This approach is 
applied to studying vibrational pure dephasing of ground state molecular iodine in a rare gas matrix. 
The Feynman-Kleinert optimized harmonic approximation for the full system density operator is 
used to sample initial conditions for the bath degrees of freedom. This extremely efficient approach is 
compared with alternative initial condition sampling techniques at low temperatures where classical 
initial condition sampling yields dephasing rates that are nearly an order of magnitude too slow 
compared with quantum initial condition sampling and experimental results. 



I. INTRODUCTION 

Detailed analysis of the vibrational spectroscopy of 
chromophores in solution can be used to investigate in- 
termolecular interactions in condensed phases. The the- 
ory of vibrational pure dephasing and its contribution 
to spectral line shapes and shifts has been worked out 
in detail in various limits. Perturbation theory (see 
references'^ and the and the literature cited therein), 
which assumes that the interaction between vibrational 
degrees of freedom and the environment is weak, gives 
an expression that enables the pure dephasing time to 
be computed from the zero frequency component of the 
time correlation function of the fluctuations in the energy 
gap between the vibrational levels. In many experiments, 
however, there may be strong initial environmental inter- 
actions associated with how the system is initially excited 
so techniques to study dephasing and dissipation beyond 
the limits of perturbation theory are important, and sev- 
eral non-perturbative ways to study vibrational dynamics 
have appeared in the literature^^ 7 ^. One very fruitful 
way to go beyond perturbation theory, for example, is 
to employ an idealized model which can be worked out 
analytically such as a two level system appropriately cou- 
pled to a harmonic bath for which the effects of environ- 
mental dephasing on lineshape and spectral shifts have 
been worked out in detaili£>£. Several general predictions 
emerge from this analysis concerning, for example, how 
lineshapes are affected by characteristics of bath spectral 
density, and temperature, etc. Experimentalists can use 
the predictions of this theory to interpret their findings 
in terms of the nature of the underlying interactions in 
condensed phase systems. The interactions between the 
vibrational coordinate and the environment can depend 
on vibrational excitation in a complex way and this can 
complicate use of the two level system theory described 
above. Also employing a harmonic bath when anhar- 
monicity may be significant could make the use of such 
a theory questionable. 

An alternative non-perturbative method that has re- 
ceived considerable attention in the quantum optics 
literatur o 9 ' 10 and recently has become a focus in molecu- 



lar science applications- is the so-called stochastic wave 
function approach 7 . With this method the evolution 
of the density matrix is assumed to take a simplified 
Bloch form parameterized by phenomenological coeffi- 
cients governing the decay of coherence and the trans- 
fer of population. Rather than propagating the n x n 
density matrix elements, a stochastic wave function rep- 
resented in terms of n relevant basis functions is evolved 
using statistical rules designed in such a way so that the 
reduced density constructed from an ensemble of evolved 
stochastic wave functions reproduces the evolution of the 
Bloch equations. Since this approach relies on a Bloch 
model form containing phenomenological parameters it 
thus provides an efficient, linear scaling approach for fit- 
ting experimental data to such models. As these models 
are often discussed in terms of gas phase collision pro- 
cesses, interpreting such fits in a meaningful way and 
extracting information about a microscopic mechanism 
of the relevant decoherence processes operating in con- 
densed phase systems is not straightforward. 

A further possible alternative approach is the use of a 
realistic microscopic model with vibrational state depen- 
dent interactions as has been developed in various con- 
texts over the last few year a 11 ' 12 ? 13 . Such an approach re- 
quires a quantum dynamical treatment of dephasing and 
microscopic simulation methods to address this problem 
have recently been develope d 14 i 15 ' 16 i 17 ' 18 i 19 ' 20 ' 21 . In this 
paper we extend these methods to include quantum ini- 
tial condition sampling that should accurately capture 
the relevant underlying spectral density of the realistic 
model system. Thus unlike the theories discussed above 
which arc built on use of different models (e.g. Debye or 
pseudolocal mode spectral densities), with the approach 
employed here we do not need to assume anything about 
the spectral density arising from the interactions. As 
such direct, model free, interpretation of the experimen- 
tal results is possible. 

The specific set of experiments which we will study 
with this alternative approach come from recent work 
from Apkarian and coworker o 22 ' 23 '^'^'^'^ 7 -'^^ 9 -'^ in 
which they use Time Resolved Coherent Anti-Stokes Ra- 
man Scattering (TR-CARS) to directly probe the de- 
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phasing of vibrationally excited I2 wave-packet compo- 
nents due to interactions of this chromophore with a 
rare gas matrix. This well studied system is chosen as 
a benchmark on which to evaluate the approach due to 
the availability of accurate interaction potential data. 
Also, as distinct from traditional frequency domain ex- 
periments in which lincshapes and frequency shifts are 
indirectly interpreted in terms of the evolution of the 
density matrix, these time domain experiments can be 
directly connected to the theory of decoherence which 
we outline below. 

The phenomenon of vibrational dephasing can be 
probed when a vibrational subsystem, described by 
Hamiltonian H s , with eigenstates given by H s \vq) — 
€ Vo \vq), is prepared in some coherent superposition state, 
IV'(O)) = J2v c va \ v o)i f° r example, in the presence of 
an environmental subsystem in state |x(0)) which is 
unaffected by this preparation. Thus the composite 
system initial state can be described by a separable 
product state 1^(0)) = \tp(0))\x(0)) ■ Important devia- 
tions from this separable product initial condition in the 
limit of strong system - environment coupling have been 
discusse d 31 ' 32 i 33 . If the quantum subsystem and the en- 
vironment are uncoupled, the initial coherent superposi- 
tion of vibrational states will be maintained as the whole 
system evolves, i.e. the amplitudes and relative phases of 
the component state contributions will be constant, and 
the composite system wave function will remain sepa- 
rable. In the more general situation, however, the dy- 
namics of the full system is governed by the coupled 
system-bath Hamiltonian, H = H s + + H s -b and the 
above initially separable wave function will evolve into an 
entangled state = Y,v °vo exp[—|#«]|x(0)vo) = 

J2v N \Xv N (t)) \vn) in which the amplitudes and rela- 
tive phases of the different vibrational state contribu- 
tions will change with time. In the last equality in 
the previous result we have projected the composite 
state onto the diabatic vibrational basis states |ujv), and 
the bath states that appear in this entanglement, la- 
belled according to this basis, are obtained as \\v N (t)) = 

( v n\ J2v Cv o ex P[~i-^]lx(0) w o)- The timescale for the 
variation in the phase of the different vibrational state 
components is known synonymously as the vibrational 
dephasing or decoherence time. The variation in the am- 
plitude, on the other hand, reflects the vibrational state 
population relaxation time. With this description, the 
bath states \xv N (t)) contain all the relative phase and 
amplitude information of the contributions from the dif- 
ferent vibrational basis states. With the view outlined 
here we are assuming that the details of the shapes of 
the exciting radiation field pulses can be neglected and 
that we can focus our attention on the evolution of the 
reduced density matrix. A complete description of the 
experiments would require considering details of the time 
dependent interactions of the full system with the radia- 
tion field and is beyond the scope of the current investi- 
gation of vibrational pure dephasing times. 

In the vibrational pure dephasing limit it is assumed 



that the composite system dynamics results in slow am- 
plitude relaxation of the chosen vibrational basis states 
on the fast timescale of the fluctuations in the phases 
of the different component basis states. Thus we sup- 
pose that contributions to final state amplitude that 
originate from initial states vq, different from vn, are 
negligible and we approximate the above expression as 
\Xv N (t)) « c VN (v N \exp[-j:Ht\x(0)v N }. This is equiva- 
lent to the vibrationally adiabatic approximation. The 
timescale for the vibrational pure dephasing process is 
governed by the distribution of fluctuations in the phase 
of the state \Xv N (t)) which is determined by two fac- 
tors: (1) the strength of the interactions between the 
vibrational subsystem and its environment which will in 
general depend on the particular vibrational state |vjv), 
and (2) the distribution of initial states of the environ- 
ment |x(0)). This paper develops a general approach 
for computing vibrational pure dephasing rates in con- 
densed phase systems which incorporates both the vary- 
ing strength of environmental interactions with quantum 
subsystem state, and the effect of quantum dispersion on 
the nature of the distribution of the initial environmental 
states. 

The quantum vibrational state dependent intermolecu- 
lar potential idea developed in the late 1990's^*^ 3 - pro- 
vides a framework for the representation that we employ 
here to account for the the variation of environmental in- 
teraction strength with vibrational state. Specifically the 
model of Martens and co-worker a 34 ' 35 i 36 ' 37 is adopted for 
this purpose. The main goal of the work here is to explore 
how thermal and quantum fluctuations of the initial envi- 
ronment influence vibrational pure dephasing in a model 
system with realistic interactions. Thus we will compare 
classical and approximate quantum descriptions of the 
thermal environmental initial distribution and explore 
their effects on dephasing dynamics. In addition a com- 
plete and general derivation of the Feynman-Kleinert- 
Wigner (FK-Wigner) approach for sampling quantum 
initial conditions is presented. 

The paper is organized as follows: First we outline 
a general approximate approach for treating the evolu- 
tion of the vibrational reduced density matrix. With this 
linearized approximation, quasi-classical trajectories are 
evolved from quantum initial conditions sampled from 
the Wigner transform of the initial density. The main 
methods section of the paper summarizes the FK-Wigner 
approach used for this sampling. Details of the complete 
derivation of this approach are given in the appendices. 
The final methods subsection makes comparisons of this 
approach with other published techniques. Finally in sec- 
tion IIIII we outline the implementation of these methods 
for application to computing vibrational pure dephas- 
ing rates and compare these data computed using differ- 
ent approximations with available experimental results. 
Concluding remarks are given in section ITVl 
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II. METHODS 
A. Vibrational Pure dephasing 

The density matrix formulation offers a convenient way 
to rewrite the wave function description outlined above 
in a useful form for making further approximations and 
interpretations. As noted earlier, the experiments of in- 
terest involve exciting the vibrational subsystem inde- 
pendently of the environment. The bath degrees of free- 
dom are thus assumed to be initially prepared in thermal 
equilibrium with the quantum subsystem. The quantum 
vibrational subsystem is then prepared by rapid pulsed 
laser excitation, for example, in a non-equilibrium co- 
herent superposition of excited vibrational states \ip{0)) 
as described in the Introduction. We assume that the 
non-equilibrium composite system is thus initially de- 
scribed by a density operator which is a product form 
p(0) = p s (0)pl(0), where pi is the bath part of the equi- 
librium density operator for which we will develop ap- 
proximations, and the non-equilibrium quantum subsys- 
tem density for the initially excited coherent superpo- 
sition state is p s (0) = |^(0))(^(0)|, which has various 
component operators, for example, c Vo c*,\vo)(v' \, that 
depend on what states of the quantum subsystem are 
coherently excited by the laser pulses. 

The coherently excited composite system will evolve 
from this factored initial state to an entangled state as 
a result of the coupled full system evolution as discussed 
above. This entanglement will thus be described by the 
full system time dependent density matrix with elements 
in the environmental coordinate and vibrational subsys- 
tem state representation given as (Rn, VN\p(t)\R' N , v' N ) 
= (R N ,v N \e- l6t / h p{0)e lf{t / h \R' N ,v' N ). The experiments 
of interest probe only the quantum subsystem states so 
we will study the reduced density matrix elements ob- 
tained by tracing over all the bath degrees of freedom 
Le - pTnv'^) = J dR N (R N ,v N \p(t)\R N ,v' N ). Suppose 
the preparation selects out the p s (0) — \vo){v \ compo- 
nent sub-system density operator initially, then the re- 
duced density operator matrix elements at time t will 
have the form 

Pltv'Jt) = JdR N J dR J dR' (R N ,v N \e~ l ^ h \Ro,v ) 
x(R \p L b \R' Q )(R' , v' 1 e tflt,h \R N ,v' N ) ( 1 ) 



and involve forward and backward propagator matrix el- 
ements as well as initial bath density operator matrix 
elements. 

Suppose the full Hamiltonian is 

H=^- + v(s) + ^-+V b (R)+^ s - b (s,R) (2) 
Ira 2M b 

which can be expressed in our quantum subsys- 
tem diabatic vibrational basis as H = P 2 /2M b + 
E Q ,/3 \a)h aP (^R)((3\ where h a p(R) = [e a + V b (R)]S aP + 
(a\$ s -b\f3)(R)- Here <& s - b is the system-bath interaction 
potential. In the case of vibrational pure dephasing we 
suppose that h a p(R) ~ for a ^ (3, i.e. the off-diagonal 
elements are small compared to the diagonal elements, so 
we can approximate the full system Hamiltonian by the 
diagonal form H d = P 2 /2M b + J2 a \a)h aa (R)(a\ for all 
important R. In this case there is no population relax- 
ation between our vibrational basis states so the propa- 
gator matrix elements appearing in Eq. ([T]) can be written 
in the composite bath subsystem path integral forms as 
follows: 

(R N ,v N \e- iAt ^\Ro,v ) 

r R(t) = R N 

= S VN . V0 V[R(t)]e^o[RW] (3) 

JR(0)=Ro 

and 

(R!o,v'o\e im/h \RN,v' N ) 

= S« s ,v> / V[R'(t)}e-^o< lR{t)] (4) 

where the forward path action, for example, is 
/"* 1 

Sv ,v [R(t)} - / dt'{-M b R 2 (t') 
Jo 1 

-[e V0 + V b (R(t')) + (v \<Z> s - b \v )(R(t'))]} (5) 

and a similar expression for the action along the back- 
ward path is obtained by modifying the vibrational state 
accordingly to v' . 

Combining these expressions, the reduced density ma- 
trix in the pure dephasing limit can be written as 



Pvo,V (t) = f dR N [ dRo f dR' (R Q \pt\R' ) [ m ^ V[R{t)] [ R(t) flN p[ jR ' W ] e H^o,.o[«W]-^,„,[i?'(t)]} (6) 

J J J JR(0)=R Q JR'(a)=R' a 

I 



In the above expression the bath evolution is still de- a computable expression for the time dependence of 
scribed at the full quantum level. To proceed to the reduced density matrix elements we follow various 
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author o i 15 ' 16 i ' 18 i 19 and combine the forward and back- 
ward propagators in Eq.([S]) and write the product in 
terms of mean, R(t) = [R(t) + R'(t)]/2, and difference, 
Z(t) = R(t) — R'(t), bath subsystem path coordinates 
(with a similar transformation for the bath momenta, 
P(t) = [P(t) + P'(t)]/2 and Y(t) = P(t) - P'(t)). Next, 
the phase of the integrand in Eq. ^ is expanded in the 
path difference variables. In condense phase systems var- 
ious arguments can be given to justify keeping only low 
order terms in the path differenc o 14 ' 20 ' 38 ' 39 ' 40 . Thus we 
proceed by truncating the expansion of the phase to lin- 
ear order in the path difference variables. With this ap- 
proximation the integrals over the initial difference coor- 
dinate Zq can be performed defining the Wigner trans- 
form of the initial density operator 

(%) w (Ro,P\) = J dZo(Ro + Y\P € b\ R o ~ ^)e- iP - z ° 

(7) 



If we discretize the paths appearing in Eq. ([6]) by inserting 
N resolutions of the identity in the bath subsystem phase 
space and take small time steps St = t/N, thus writing 
the discrete path variables as Rk = R(kSt), for example, 
we find that all integrals over the difference coordinates, 
Zk, and difference momenta, Y k for < k < N that 
appear in the discretized result can also be performed 
in the linearized approximation since they are integral 
representations of <5- function s 20 ' 41 . Thus, the linearized 
approximation for the evolution of the reduced density 
operator becomes 



Pv ,v' (t) 



N 



dPk 



dR YYdR k —{pl) w {R 0l P 1 )e h( "° V e 



k=l 



N-l 

fe=i 



/ p k+1 -p k 

V St 



pvo,v 



N 



fe=i 



Pk Rk — Rk-1 

M " St 



(8) 



where 

K 0,< = -\ (VhA«(4) + Vx k h v , K (R k )} (9) 
and 

A1V„< = (v Q \<l>s- b \vo) - (v'o\^s- b \v'o) (10) 

This result indicates that the time evolution of the den- 
sity operator matrix elements can be approximated by 
first sampling the initial environment phase space from 
the Wigner transform of the bath part of the thermal 
equilibrium density of the full system. In our low temper- 
ature calculations we approximate this equilibrium den- 
sity by assuming that the environment experiences in- 
teractions with only the ground state of the vibrational 
subsystem. At higher temperatures excited vibrational 
states should be included in this initial sampling. Next, 
the product of 5-functions in Eq.® indicates that with 
in the linearized approximation the time dependence of 
the density matrix elements can be obtained by evolving 
classical trajectories with the sampled initial conditions 
and subject to the mean of the forces arising from the 
quantum states involved in the prepared superposition 
state. Finally Eq.© gives that the contributions from 
each sampled trajectory to the density are obtained by 
adding coherently phase factors computed along these 



trajectories. An approach based on the same approxi- 
mations outlined above has been used in other wor k 20 i 41 
to compute various quantum time correlation functions 
exploring electron transport, and vibrational energy re- 
laxation, for example. In fact, approaches like this which 
employ the Wigner transforms 2 - of the equilibrium den- 
sity as the distribution of initial conditions and evolve the 
classical dynamics with a mean Hamiltonian have a long 
history in computing spectroscopic correlation functions, 
for example ) 43 ' 44 . The work of Egorov et a?4& provides 
an important comparison of this type of approach with 
alternative classical and mixed quantum-classical meth- 
ods in the context of computing model condensed phase 
vibronic spectra. 

The rest of the paper describes the computational ap- 
proach we adopt for sampling the Wigner density and 
the mean Hamiltionian dynamics calculations that we 
have conducted to explore the characteristics of the vi- 
brational pure dephasing process in low temperatures 
crystals. In Fig. [1] we present our dephasing results 
for I2 in solid krypton and compare with experiments 
in order to demonstrate the sensitivity of these type of 
calculations to the potential field used in the underlying 
dynamics. The calculation results presented in this fig- 
ure both use the FK- Wigner initial condition sampling 
approach detailed in the next section. The different sets 
of calculated dephasing rates presented here, however, 



5 



0.12 

> °- 1 
o 0.08 

H — ' 

CO 

I 5 0.06 

c/5 

CO 
JZ 
CL 
CD 

Q 



0.04 



0.02 



FK-Wigner IC ASF T=2.6K 
FK-Wigner IC GSF T=2.6K 
Experiment T=2.6K-' 



2 4 6 8 10 12 14 
l 2 vibrational state number v 



16 



FIG. 1: Comparison of experimental and computed dephasing 
rates versus vibrational state for 12 in solid Krypton at low 
temperatures. Solid curve gives results obtained using trajec- 
tories propagated with the Average vibrational State Force 
(ASF), dashed curve gives results obtained with vibrational 
Ground State Force (GSF) (see text for discussion of these 
different approaches), and dotted curve displays the experi- 
mental results. 



were obtained using different ways of propagating the 
classical trajectories. The solid curve presents results 
obtained using the average state force (ASF) to drive 
the environmental dynamics as suggested by Eq.© and 
discussed above. An alternative approach is to evolve 
the environmental dynamics using forces arising from the 
ground vibrational state potential onl y 43 ' 45 . The curve 
labeled GSF (Ground vibrational State Force) in Fig. [1] 
is obtained in this way. We see that both sets of cal- 
culated dephasing rates are generally faster than the ex- 
perimentally observed results. These discrepancies with 
experiment probably reflect inaccuracy in the model in- 
teractions employed in these calculations which we out- 
line later. However, we generally find that the dephas- 
ing rates computed using the average surface force are in 
better agreement with experiment than results computed 
using the ground vibrational state force. The two calcu- 
lated dephasing rate curves should coalesce at low vi- 
brational state number where the average potential from 
states and v is essentially that of the ground state. 
As we increase v the average force experienced by the 
environment = [V R h . (R) + \7 R h^ u (R)]/2 in our 
ASF calculations becomes more and more different from 
the ground state force F 0,0 = V^/iq,o(-R) employed in 
the GSF calculations and the deviation between the two 
curves is apparent. As discussed above, according to Eq. 
([8]), the dephasing rates are obtained by averaging the 
phase factor in the energy gap between the two relevant 
vibrational states. If the forces governing the fluctua- 
tions in the energy gap arise from both of the states we 



expect that each state energy will fluctuate in a similar 
way leading to slower dephasing. If, on the other hand, 
the solvent responds to only the ground state force, as in 
the GSF calculations, the excited state energy may fluc- 
tuated more strongly as it does not play a role in influ- 
encing the environmental dynamics. This explains why 
the dephasing rates computed with the average force dy- 
namics are slower than those obtained with the ground 
state force. 

The main goal of our work is to study the effect of 
the distribution of environmental states on vibrational 
pure dephasing, making particular contact with low tem- 
perature experiments. Under these conditions quantum 
dispersion and tunneling of the environmental degrees of 
freedom may play a significant role. Thus we will em- 
ploy various approximations to the Wigner transform of 
the initial density and compute their effect of the lin- 
earized approximation to the dephasing dynamics given 
in Eq.©. General methods for computing the Wigner 
transform of the Boltzmann operator are as yet unavail- 
able. In our work we thus employ an approximation 
to this operator. First we assume that the tempera- 
ture is sufficiently high, and that particle masses are 
large, so that identical particle statistics complications 
can be ignored. Next, as detailed below, we employ an 
approach pioneered by Feynman and Kleinert that ap- 
proximates the many body Boltzmann operator assum- 
ing a locally quadratic form for the interactions. This 
approximate approach has recently been implemented 
in various condensed phase applications studying quan- 
tum effects on transport and spectroscopy by Poulscn 
and co- worker s 14 i 46 1 47 . Our study will compare vibra- 
tional dephasing results obtained using this approximate 
Boltzmann operator which can incorporate some quan- 
tum initial distribution effects, with classical initial con- 
dition sampling techniques which ignore quantum dis- 
persion and tunneling. Our results from all these various 
calculations will be compared with available experimen- 
tal results enabling a detailed understanding of the reli- 
ability of the different approximations underlying these 
calculations. 



Sampling the quantum initial density for the 
environment 



In this section we review the method developed 
by Poulsen and co-workers 1 ^^^^^^^^^ that 
adapts the Feynman-Kleinert variational approach for 
obtaining an optimal local harmonic approximation to 
the Boltzmann operator which can be easily Wigner 
transformed. In our application to vibrational dephas- 
ing we suppose that the quantum vibration is prepared 
initially in some superposition state while the environ- 
mental degrees of freedom are unaffected by the vibra- 
tional excitation and are initially in thermal equilibrium 
with the vibrational degrees of freedom. As outlined in 
the previous section this quantity plays a central role as 
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the phase space distribution function that must sampled 
to provide trajectory initial conditions in linearized ap- 
proximations for the density matrix dynamics. We will 
explore the effects of the quantum nature of this initial 
equilibrium distribution on vibrational pure dephasing of 
the excited coherence. The development below in terms 
of the general Boltzmann operator is easily adapted to 
the case of the environmental degrees of freedom being 
in equilibrium with the quantum subsystem by supposing 
that the vibrational degree of freedom acts to provide an 
external force on the environment. As mentioned above 
we suppose that the experiments of interest can be inter- 
preted as though the environment was initially in equi- 
librium with the ground state of the quantum subsystem 
so the Hamiltonian in the Boltzmann operator below de- 
scribes interactions of the environment with the external 
field due to the ground state quantum subsystem. 

The principle assumption underlying the Poulsen et 
a^i 14 i 46 i 47 approach to generating an approximate Wigner 
transform of the Boltzmann operator, p — exp[— /3-ff], is 
that the off-diagonal elements of the thermal density ma- 
trix pp(R, R') that enter the Wigner transform expression 
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pJ(R,P) = J dAR(R + AR/2\ Pl3 \R~ AR/2)e-i 

(11) 

can be just as well represented as the diagonal elements 
that appear in the trace expression for the partition func- 
tion. Accepting this proposal the approach proceeds as 
follows: First the Feynman-Kleinert method for comput- 
ing the partition function is employed to obtain a local 
harmonic approximation to the Euclidian action that is 
variationally optimized for computing the trace of the 
density matrix. Next we assume that the off-diagonal 
elements can be approximated using the optimized lo- 
cal harmonic approximation and this form is employed 
to compute the Wigner transformation. With this ap- 
proach a local harmonic approximation is never directly 
made to the potential, rather an optimal local harmonic 
approximation to the Euclidean time action is found vari- 
ationally using the partition function which involves in- 
tegration over all space and thus contains global infor- 
mation. This local harmonic fitting procedure, however, 
is not influenced by any off-diagonal information. So it 
is unclear why it should provide a reliable approxima- 
tion for generating the Wigner phase space density. One 
of our goals is to explore this question. In section III CI 
we will directly compare density matrix elements for a 
simple model computed with in this approximation with 
exact results to explore the reliability of the off-diagonal 
density obtained with this and other approximate meth- 
ods. 

As outlined above the derivation of the approach 
starts by following Feynman and Klcinert and consid- 
ering the partition function, Z , at inverse temperature, 
(3 = 1/ksT, for the environmental variables, R, writ- 
ten as a path integral over cyclic paths (R(t) : R(0) = 



R((3h)) i.e. 



Z = Iv[R(r)]e 



-S[R(r)]/h 



(12) 



Where, §T>[R(t)], means integrate over all points in 
the path, i?(r), including the common starting and fin- 
ishing points, and the Euclidian action is S[R(t)] = 

^ h dT{\R T { T )M.R{T) + V{R{T))}. Here, M, is the mass 
tensor, and V(R), is the inter-particle interaction poten- 
tial including the effects of external fields (such as, for ex- 
ample, coming from interactions with the vibrator in its 
ground state). In mass weighted cartesian coordinates, 
q(r) = M 1/2 R(t), we can write 

Z=l V[q(r)}e-^^ h dr[^(r) + v {q{ r))] (13) 



Since the paths of interest for Z are real and periodic 
on < r < (3h their Euclidean time dependence can be 
written in terms of a Fourier series i. e. 



q(r) = q c + > </„< 



]>>„ e *°" T + q* n e- iQ ^\ (14) 



where f2„ = 27tn/f3h are the Matsubara frequen- 
cies, the real zero frequency Fourier coefficient, q c = 
~BK f() h ^ r< ?( T )' i s l ne P a th centroid, and all other Fourier 
coefficients, q n , are in general complex. 

If the configuration of the system is of dimension d then 
the cyclic path space integration can be mapped onto an 
integration over Fourier coefficient vectors of length d 
according to£££&£££& 



V[q{r)} 



dq c 



{2nh 2 f3) d /< 



n j dRcq n J dlmq n 



and the partition function can be shown to be 



dq c 



n 



/ dReq n J dlmq n 



xe 



E£U a* \q„ | 2 e - ±Jt drV(q(r)) 



(16) 



Following Feynman and Kleiner t 55 i 56 the effective classi- 
cal potential, W(q c ), is defined as 



-pw(q c 



n 



J dReq n j dlmq n 



X£r /3£~ =1 ^M 2 e ~i JT drV( q (r)) ^ 



and the quantum partition function is cast in the sugges- 
tive classical form 



Z = 



dq c 



(2irh 2 P) d / 2 



(18) 



identifying motion of the path centroid over the ef- 
fective potential as a classical-like way to compute 
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quantum properties using classical statistical mechanics 
ideas^£££2. Since (5n 2 n = Air 2 n 2 /h 2 f3, at high tempera- 
ture (lim/3 — > 0), the Gaussian factors in the integrand 
of Ea. fTT]) , for example, will be dominated by small val- 
ues of \q n \> so under these conditions q(r) will fluctu- 
ate only slightly from q c and, as Feynman and Kleinert 
suggested^, the effective potential might be computed 
using a local harmonic approximation to the full potential 
in the region around each path centroid location. They 
used the Gibbs-Bogoliobov-Feynman variational princi- 
ple to determine the optimal local harmonic approxima- 
tion at each centroid position. 

This approach starts from the exact expression for 
the partition function obtained by multiplying and 
dividing by some trial partition function Z tr iai = 
§ T>[q(T)\e~~n Stl ' ial \- q ( T ^ associated with some, yet to be 
specified, locally harmonic trial action, S tr i a i, thus 

§V[q(T)]e-h{ S [9(-r)]-^rial[g{r)]} e -^S Mal [q(r)] 
Z = ; : — — — , , Ztriai 



e -^{S[g(r)]-S trIa! [g(r)]}\ ^ 



/ 



trial 



trial 



(19) 



Due to the concavity of the exponential function, (e~?) > 
e~^\ so we replace the average of the exponential by 
the exponential of a more simply computed average and 
obtain the following variational result 



Z > e -\{S[q(r)\-S t „ al [q(r)\) t 



J trial 



(20) 



If we choose St r ial[q{r)} = J dT{\q 2 {r) + V tr i a i(q(T))} 
to be the Euclidian time action associated with motion 
in some trial potential Vtriai-, then 

(S[q(r)] - S t rial[q(l~)]) trial 

= /J m dr[V{q{r)) - V trial (q(r))]\ (21) 

\ / trial 




dr'f(q(r')) 



trial 



Ztriai JQ 



dr' <b V[q(T)]f(q(T'))e-i St "^ T \k2) 



and the value of the cyclic path integral should be inde- 
pendent of the value of the time at which the function / 
is evaluated since all r' should be equivalent under the 
trace we can write^, for example, 




dr'f(q(r')) 



ph(f(q(0))) tnal (23) 



trial 



To proceed we follow Feynman and Kleinert and suppose 
that the trial action is quadratic in displacements about 



the path centroid location thus 
f ph 1 

Strial[q{r)] = J dT-[q{T) T q{T) 

+ (q(r) - q c ) T K(q c )(q(r) - q c ) + L(q c )} 

(24) 

The term linear in displacement, (q(r) — q c ), multiplying 
a derivative vector evaluated at q c vanishes identically 
when integrated over r due to the cyclic nature of the 
paths on the interval < r < 0H fEq. (fT4"]) ), so with this 
quadratic form only the offset term, L(q c ), and the local 
curvature matrix, K(g c ), need to be determined varia- 
tionally. 

Suppose the matrix, U(g c ), diagonalizes the curvature 
matrix giving a diagonal matrix of frequencies, u)(q c ), 
according to U T KU = ui 2 . Transforming to normal 
mode vectors, ?/(r) = U T (j(t) 



1 we find that the trial partition function is ob- 



tained as 



atrial 



drj c 



n 



Dn functio 
J dRer] n J dlmrj^ 



The Gaussian integrals over the components of the r\ n 
can be performed analytically provided the frequencies 
u)j have appropriate values (see detailed discussion below 
on the restrictions that this places on the matrix K) and, 
using the fact that 



sinh : 



x 



= l[[l + X 2 /(n 2 n 2 )} 



(26) 



the trial partition function becomes 

d 



Ztriai — 



dR r . 



(2TTh 2 P) d / 2 



n 



11 8wh./3hu)j(Rc)/2 

(27) 

so in analogy to Eq. (|17[) we can write 



WtrURc) - -^ln ( vSinh/3 ^ (i?c)/2 



L(R C ) 



The other quantities that need to be computed to ap- 
ply the variational result in Eq. ([20|) according to Eqs. ([2Tj) 
- (123)) , are averages of the potentials over the trial den- 
sity. In Appendix [Al we outline the computation of these 
quantities obtaining the following general result written 
in terms of different Gaussian smeared potential forms: 



(V(R)) 



1 



trial 



dR r 



ZtrialJ (2lTh 2 (3) d / 2 
d 

X 



n JMft. *w (») 



sxah/3tTWj(R c )/2 
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where Va{R c ) is the Gaussian smeared full potential 



Va(Rc) 



dR 



|27rA|V2 



■±{R-R C ) T A- 1 (R C ){R-R C ) 



V(R) 
(30) 

which can also be computed in Fourier space according 
to Eq. (|A7[) . The Gaussian smeared local harmonic ap- 
proximate trial potential is 



V A nal (R c ) 



dR 



,-±(R-R c ) T A- 1 (R c )(R-R c ) 



, |2ttA|V2 

i? c ) T M 1 / 2 X(i? c )M 1 / 2 (i? - R c ) + L(R C ) 
WM^K^M^A^+LiRj 



(31) 



which can be used in Eq. (f2"5)) to compute 
(Vtrial(R(fy)) trial' bi these results the smearing 
width matrix A is 



A = M- 1/2 UAU T M- 1/2 



with 



{ m 2 n+^(Rc)} 



(32) 



(33) 



as obtained in Appendix O (See Eq. (|A"5l) and Eq. (fA"6j) ). 

Thus we find the trial averaged action difference ap- 
pearing in the variational expression, Eq. (|20p . can be 
written as 



(S-S t rial) trial trial 

d 

3 



1 



dRr. 



PURc) J~J Phuij{R c )/2 



x 



ZtrialJ {2irh 2 f3) d / 2 lJ { smh{3huj J (R c )/2 
ph[V A (R c ) - ~ Ml /2 K t3 A tJ M^ /2 - L(R C )} (34) 



Using the above result we follow Feynman and 
Kleinert5£ and optimize the right hand side of Eq. (f2"0l . 
/ = exp[-i(5 - Strial)trial] z trial,'wtih respect to the 
variational parameter functions K(i? c ) and L(R C ). The 
details of this variational calculation are given in Ap- 
pendix [5] and it gives the following relationships between 
the different optimized parameters: 

L(R C ) = Va(R c ) -\Y, Ml /2 K tJ (R c )A tJ (R c )My 2 

(35) 



Here the optimal curvature matrix is found from the 
gaussian smeared result 



K<KJ = /i^ D(ii » e 

Where D(i?) is the mass weighted Hessian 



\(R-R C ) T A- 1 (R C )(R-R C ) 

(36) 



1/2 



(37) 



As outlined above the approach requires the computa- 
tion of Gaussian smeared averages of various functions. 
The evaluation of such multidimensional integrals is in 
general complicated and could be performed with a MC 
procedure. In the application described here, however, 
the interaction potential functions being smeared are pair 
decomposable and, as outlined in Appendix |C| can be fit 
to Gaussian forms and the smearing integrals performed 
analytically. 

The structure of the above equations suggests the fol- 
lowing iterative procedure for computing the smearing 
matrix A: (1) Choose an initial guess for the A( ) matrix. 
We use the previously converged smearing matrix for the 
last accepted centroid configuration. (2) Compute the 
full K( 0) matrix at the new centroid configuration using 
Eq. (|3"6")) (or Eq. (|C4[0 . This can be done efficiently em- 
ploying the symmetry of the local curvature matrix. (3) 
Diagonalizc K^ ) to obtain the local harmonic frequen- 
cies u/°)(i? c ) and the eigenvector matrix U^ -* as defined 
above Eq.([n5]). Finally, (4), use these in Eqs.([32]) (or 
(|B16p ) and (|32[) to obtain a new estimate of the smear- 
ing matrix, A*- 1 -* , then return to step (2) and iterate this 
process till convergence. 

Using the converged set of parameters, we compute 
the various Gaussian smeared quantities (using Eq. (f30"|) 
for example) and obtain the variationally optimized value 
for the energy offset parameter L(R C ) as in Eq. ([33|) at the 
new centroid configuration. 

The final steps in the approach involve writing the den- 
sity operator approximately as the integral over centroid 
locations of the variationally optimized local harmonic 
form around each centroid position and then analytically 
Wigner transforming this locally harmonic result. De- 
tail of these calculations are given in Appendix [D] The 
approximate density has the form 
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x 11^ { ^ - fn ™^/2M - ,D 2 } /^/S (38) 

Here we define the quantity aj = coth fihwi /2 — 2/f3hiUi and Wigner transformation yields 

x expl-^fO, - mp| _ t.nh(/i/ to ,/2) 



The integral in Eq. (|39p is easily computed by impor- 
tance sampled Monte Carlo. Thus a Cartesian centroid 
configuration, R c , is sampled. The procedure outlined 
above Eq. ([55]) is employed to iterate the width, A^ n \ 
and curvature, K>"), matrices to convergence. As in- 
dicated, this involves diagonalizing KV 1 ' at each itera- 
tion to obtain the local harmonic frequencies and nor- 
mal modes. The converged value of W tr iai computed 
from Eq. (j28f is used to accept or reject the sampled 
centroid configuration. Next, to sample the (Q,Pq) 
phase space distribution associated with the approxi- 
mate Wigner transformed Boltzmann operator, a set of 
k = 1,M normal mode displacements with components 
£i = Qi ~ 7 ld are sampled from the component Gaus- 
sian distributions with variances uqi = (aih/2u;i) 1 / 2 
and the set of M points in Cartesian space B> k ' = 
R c + M _1 / 2 U(i? c )^ fc ^ is generated to provide a sam- 
pling of initial configurations around the point R c . Simi- 

(k) 

larly, Gaussian random numbers Q are sampled with 
variance opi = (huji/2 tanh(/3?icji/2)) 1 / 2 to provide a 
set of normal mode momentum vectors Pq^ which are 
transformed to cartesian initial momenta according to 
P (fe) = M 1 / 2 U(i? c )C ( ' c) . The numerical factor (. . . )*/ 2 in 
Eq. (|39[) provides an automatic centroid configuration de- 
pendent normalization for the sampling of these Gaussian 
distributions, and so each sample phase space configura- 
tion generated in this way carries unit weight. 



C. Comparison of Initial Condition Sampling 
Methods for Model System 

To test the reliability of our implementation of the 
Feynman-Kleinert approach for sampling the Wigner dis- 
tribution outlined in the previous subsection we have ap- 
plied it to several simple exactly solvable one dimensional 
model systems. In this section we present these results 
and compare with exact calculations and with another 
approximate approach due to Shi and Geva (SG) 15 i 16 ' 61 . 
These authors proceed by multiplying and dividing the 



Wigner transform expression of the Boltzmann operator 
by the diagonal elements of the thermal density operator, 
thus 



(40) 



(e-^ B ) w (Q,P) = (Q\e-P B \Q) I dze' lPz ' h 

{Q-z/2\e~P H \Q + z/2) 
(Q\e-P H \Q) 



They make the Local Harmonic Approximation (LHA) to 
the potential, expanding to quadratic order about the po- 
sition Q. Note this approximation is made only in the ra- 
tio of off-diagonal to diagonal density matrix elements in 
the integrand, however, the full anharmonic dependence 
of prefactor diagonal element is included through path 
integral calculations. With this local harmonic form, the 
Gaussian integrals can be performed analytically yielding 
the following result for the case of one dimension: 



(e-? H ) s w G (Q,P) = (Q\e- pH \Q)( 
x exp[ 



Airh 



Mco(Q)x(Q) 
P 2 



hMu{Q) X {Qy 



)l/2 

(41) 



Where x(Q) depends on the local curvature of the po- 
tential about the point Q and has the form x(Q) = 
coth with w(Q) = [d 2 V(Q)/dQ 2 /M] 1 / 2 . 

For the purpose of comparison with exact results we 
computed the density matrix for our ID model using the 
numerical matrix multiplication (NMM) approac h 62 ' 63 . 
It is most convenient to make direct comparison with the 
full density matrix. The approximate Wigner densities 
can be inverse Wigner transformed to give approximate 
coordinate space density matrices using to the following 
result: 

{Q\e- 0H \Q') = {2-kK)- 1 [ dP 



xe-* p < Q - Q 'V h (e-P H ) w (-(Q + Q'),P) (42) 
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Thus the SG approximation to the density matrix be- 
comes 



-Ph^sg _ ,n\„-PH\n\ ,W- M ^)*(Q) (Q _ Q , )2] satisfy > -(2tt//3») 



Eg. ([32|) the eigenvalues are proportional to An as given 
in Eq.lO (or Eq.((BIBJ)). The right hand side of this 
expression is positive provided all the eigenvalues of K 



(Q|e-HQ'r G = (Q\e- pM \Q) exp[ 



4fi. 



K* K i.e. some negative 



(43) 

with Q = (Q + Q')/2. The corresponding FK approxi- 
mation to the density matrix takes the form: 



(Q\e- pH \Q') 



2M 2 lu(Q c 



\FK 



dQ c e 



(3ha 



1/2 



exp (Q - QcY 

ha 



x eM _M^Q c )\«m+2,^ m]{Q Q , m 



The FK form in Eq. (|44]) and the SG approximate den- 
sity matrix of Eq. (|43p share some similar features. Set- 
ting Q — Q' in Eq. (|4"4")l gives an expression for the diag- 
onal elements (Q\e~P H \Q) FK . If we evaluate the gaus- 
sian in (Q — Q') of Eq. (|4"4")) at the maximum of the first 
gaussian in this equation i.e. Q c = Q in an effort to ap- 
proximate the Q c integration, we recover a form similar 
to the SG approximation of Eq. (j43|) . In general however 
the full integral over the centroid position must be per- 
formed so we expect differences between these results. 
The most significant differences must be due to the fact 
that the SG form employs a single local harmonic fre- 
quency u>(Q) to compute the off-diagonal elements. The 
FK expression, on the other hand, combines results from 
a range of smeared frequencies that are not obtained by 
local harmonic approximation to the bare potential. 

These differences are apparent in the results presented 
in Figures. [2 and [3l In these calculations we explored 
the approximate density matrixes for an asymmetric dou- 
ble well potential (the Veneziano potential) of the form 
V{Q) = \E C (1 + Q 2 ) 2 (l - mxf with the following pa- 



rameter values in atomic units: E c — 1 x 10 



0.2 



and M = 1600. Results are presented with a temperature 
corresponding T=50K. 

From Fig. [2] it is clear that under these conditions the 
FK approximation reproduces the density matrix very 
accurately across the entire region, even when barrier 
tunneling is important. In Fig. [3j on the other hand, we 
see that the SG approach gives undefined results for the 
density matrix in the tunneling region. Moreover, the 
density matrix difference results presented in the right 
hand panel show that even though the SG density matrix 
is exact along the diagonal (when it can be defined) the 
off diagonal elements differ considerably from the exact 
results compared to the significantly smaller differences 
observed with the FK approximation. 

The general condition for convergence of the various 
multidimensional Gaussian integrals like Eq. (|A7[) in k- 
space, or analogously Eq. (f30|) in R-space is that the 
gaussian width matrix, A, appearing in these expres- 
sions should be positive definite, i.e., all its eigenval- 
ues should be positive. From the definition of A in 



frequencies can be tolerated with in the FK approach 
and this is an important factor in reproducing the den- 
sity matrix in the tunneling region. Analysis of the SG 
density matrix expression in Eq. (|43p reveals that it too 
can tolerate negative curvature regions. It is found that 
uj 2 > — (7r//3ft) 2 = —K^ 15 . Thus the maximum nega- 
tive curvatures K^ K that the FK approach can tolerate 
are 4 times those of the SG approach. Moreover, due to 
the use of the smeared potential in the FK approach as 
distinct from the bare potential in the SG formulation, 
curvatures with the FK approach are generally signifi- 
cantly smaller in absolute value than those of the SG 
approach. Thus the FK approach is expected to have a 
considerably wider range of applicability. 



III. RESULTS 

A. Outline of Experiments 

The TR-CARS experiments of Apkarian and co- 
worke ra 22 i 23 i 24 i 25 i 26 . i 2Z i 2£ i 2£ i 2£ involve exciting controllable 
coherent superpositions of I2 vibrational states in rare gas 
matrices. After the excitation pulses that prepare the ini- 
tial vibrational superposition on a chosen electronic state, 
probing pulses project the evolving packet onto other 
electronic states and the time dependence of the emis- 
sion from these states gives a signal that can be related 
to the evolution of the initial coherence. Apkarian and 
co-workers use a model of these experiments to extract 
dephasing rates of the different vibrational superposition 
states they prepare. Their work exploring these dephas- 
ing rates for vibrational superpositions prepared in the 
ground X electronic state is the focus of our studies here. 
In our work we assume that the pure dephasing rate can 
be obtained from the long time exponential decay rate of 
the off-diagonal elements of the density matrix in an ap- 
propriately chosen vibrational representation. Martens 
and co-worker a 34 i 35 i 36 i 37 have used their semi-classical Li- 
ouville dynamics approach (which is formally equivalent 
to the linearized dynamics employed here) to study this 
problem and have presented a useful model of the rel- 
evant vibrational state dependent interactions. We fol- 
low these workers and assume that the eigenstates of a 
Morse oscillator whose parameters are fit to give solu- 
tion phase experimental vibrational data for the ground 
X electronic state provides such an appropriate vibra- 
tional representation. Using this model of the interac- 
tions we explore the linearized approach for computing 
the dynamics of the vibrational density matrix with dif- 
ferent distributions of initial conditions. The main dif- 
ference between our calculations and this previous work 
is the use of the FK-Wigner quantum initial condition 
sampling approach developed in the previous section. 
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FIG. 3: Left panel: Asymmetric double well density matrix, green surface is exact numerical calculation and red surface are 
results from SG approximation. Right panel: Red surface gives difference between exact density matrix and FK approximation. 
Green surface presents differences between exact density matrix and SG approximation. 



B. Computational Model 

The computational model employed in our studies of 
vibrational dephasing is closely related to that of Martens 
and co-worker a 34 i 35 i 36 i^ and Meier and Beswic k 65 ' 66 . 
Thus we write the Hamiltonian as H = H s + Hb + H s -b 
where 



(45) 



states of system Hamiltonian, i.e. H s \v) = e v \v). We 
employ the approach of Martens and co-worker o 35 ' 36 i 37 
and approximate the density associated with these vi- 
brational basis states using two weighted (S-functions. 
Thus we write the vibrational wave functions density as 
|(r|?;)| 2 = X)fc=i c V&{r— r™) and the 5-function positions 
r™ and weights c™, which are state dependent, are fit to 
gas phase vibrational bond length moments^. Here we 
explore vibrational pure dephasing so h a p for a/^ are 
assumed to be small and the full Hamiltonian becomes 



H s -b — 



2/xf 



+ V s - b {r,0,(j),f C om,Q) (47) 



0.2 



Here \i is the reduced mass of the diatomic whose posi- 
tion is specified by f, 9, 0, f com and Q describes the con- 
figuration of the bath. In our calculations we represent 
the quantum vibration in terms of a basis set of eigen 



H = 



2m 2M 



p2 



l,i,L,Q)(a\ (48) 



where the diagonal Hamiltonian matrix elements are ob- 
tained by summing interactions of the different bond 
length molecular representations with the environment 
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as 



h aa (f com ,d,<j>,L,Q) = e v + V b (Q) 

,,'2 




sin 2 6 J 2/*E*cT»T 

k 



(49) ~ 



> 
o 
Q. 



In our calculations we employed the Morse potential 
model of Martens and coworker o 35 ^ 36 ' 37 to describe the I2 
vibrator and the interactions between the iodine atoms 
and the krypton particles as well as the solvent-solvent 
interactions are modeled using the Lennard- Jones poten- 
tials of these workers. 

This Hamiltonian model is incorporated in Eq.((8]) 
which gives that the reduced vibrational density matrix 
elements can be computed by averaging the basis state 
energy gap phase factor along a classical trajectory in 
which the environment moves on the mean potential sur- 
face produced by the two states involved in the coher- 
ence. In our studies the initial conditions for the envi- 
ronmental variables are sampled from either the classi- 
cal distribution (as used in Marten's earlier work) or the 
FK approximation to the Wigner transform of the bath 
density matrix. In our implementation of this quantum 
initial condition sampling the rotor is fixed at the mini- 
mum energy orientation and center of mass position for 
the ground vibrational state and optimized solvent ge- 
ometry. The solvent particle positions and momenta are 
then sampled from the approximate Wigner transform 
in the presence of the fixed rotor. This approach thus 
neglects the quantum rotational dispersion. The rota- 
tional dynamics, however is incorporated, as the system 
evolves on the mean potential surface using rigid rotor 
ME)S with the weighted moment of inertia appearing in 
the second last term in Eq . (|4"9")l . 

In our studies the iodine molecule replaces two krypton 
atoms in a double substitutional site in a 108 particle 
FCC lattice. 

In implementing the FK sampling approach of Eq. |59"1) 
we perform regular Metropolis MC sampling of the cen- 
troid position and accept configurations based on the 
change in Wtrial- The solution of the variational equa- 
tions was generally found to converge in a very small 
number of iterations, typically two or three. For each 
sampled centroid position we use gaussian sampling to 
generate 5 phase space points according to the FK ap- 
proximate Wigner distribution which serve as trajectory 
initial conditions. All the results presented below have 
been averaged over 5000 independent trajectories. 



C. Dephasing Time Results 

In Fig. 0] we present results showing the time depen- 
dence of the off-diagonal elements (ground and excited 
state) of the reduced vibrational density matrix. As ex- 




FIG. 4: Decay of various ground - excited state elements of the 
reduced vibrational density matrix of I2 in for solid krypton at 
T=2.6K obtained with FK Wigner initial condition sampling. 



pected the coherence decays extremely slowly for exci- 
tation of low energy superposition states. However, due 
to the disparate nature of the interactions between the 
ground state molecule with its environment and a highly 
excited molecule and its environment, the large fluctua- 
tions in energy gap between these states, in this situation, 
cause rapid decay of these off-diagonal density matrix el- 
ements. 

We have extracted the dephasing rates from our lin- 
earized dynamics calculations of the reduced density ma- 
trix by fitting the long time behavior of the data in Fig. [4] 
to an exponential decay form. The computed dephasing 
rates, for classical and FK- Wigner initial condition sam- 
pling calculations are compared with the experimental 
results of Apkarian and co-worker a 28 i 29 in Fig. [SJ Here 
it is clear that at the higher temperatures, T=20K and 
32K, (center and right hand panels) the nature of the ini- 
tial condition sampling has little effect on the calculated 
dephasing rate and that these calculation results agree 
well with the experimental results. The fact that the re- 
sults from FK- Wigner and classical initial condition sam- 
pling agree reasonably well at high temperatures suggests 
that our implementation of the FK- Wigner sampling ap- 
proach is reliable. 

At low temperatures, for example T=2.6K, however, 
(see left panel in Fig. [5]) the approach used to sample 
initial conditions has a significant influence on the com- 
puted dephasing rates. The rates obtained with classi- 
cal initial conditions are almost and order of magnitude 
smaller than the experimental dephasing rates observed 
at this low temperature. Dephasing rates obtained with 
the FK- Wigner sampling approach, however, are on the 
same order as experimental results and show a very sim- 
ilar increase in dephasing rate with vibrational state. As 
discussed in detail by Apkarian and co-worker o 28 ' 29 , these 
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low temperature experimental results presented in Fig. 
[5] are influenced by strain induced inhomogeneity of the 
trapping sites. They suggest that their results at the low 
quantum numbers are thus sensitive to sample prepara- 
tion at these low temperatures. The fact that the slope 
of the experimental curve agrees well with our computed 
dephasing rates when FK-Wigner initial condition sam- 
pling is employed suggests that this approach is a reason- 
ably accurate way of incorporating quantum dispersion 
effects in linearized path integral dynamics calculations. 
While showing good qualitative behavior, generally the 
dephasing rates obtained with FK-Wigner sampling are 
about a factor of 2 too fast. This discrepancy could arise 
from two possible sources: First the use of the FK varia- 
tional approach for approximating the Wigner transform 
means that the trial harmonic free energy is always larger 
than the true free energy. Thus the approach overesti- 
mates the effects of quantum fluctuations and dispersion 
which is essentially equivalent to a higher effective tem- 
perature and could result in too rapid dephasing. The 
second source of deviation is the semiclassical nature of 
the linearized dynamics in which classical trajectories are 
evolved over the mean surface from the approximate ini- 
tial condition distribution. 

In the next section we explore the reasons behind the 
difference between classical and FK-Wigner dephasing 
rates by comparing structural differences in the simu- 
lated environments. 



D. Equilibrium Structure 

In Fig. [5] we compare pair distributions computed 
with FK-Wigner equilibrium distribution sampling using 
Eq. ([55)) and results from classical Boltzmann distribution 
sampling in a pure Krypton crystal. At T=32K the clas- 
sical and quantum pair distributions agree very closely 
suggesting that the temperature is sufficiently high under 
these conditions that quantum effects are unimportant. 
These distributions provide initial conditions for the dy- 
namics underlying the dephasing rate calculations. Thus 
the similarities in structure manifest themselves in simi- 
lar dephasing rates as presented in the right panel of Fig. 
[5j At T=2.6K, however, the classical pair distribution 
function is strongly peaked compared to the FK-Wigner 
result (see right panel in Fig. [S]) so the distribution of 
trajectory energy gap fluctuations resulting from classi- 
cal sampling will be narrow, leading to slow dephasing 
compared to the FK-Wigner result as discussed in con- 
nection with Fig. [5j 

In Fig. [7] we see that the quantum structure is fairly 
insensitive to a more than ten fold increase in temper- 
ature (T=2.6K to T=32K). The width of the pair dis- 
tribution function peaks in this figure increase monoton- 
ically by about 20% over this temperature range. Fig- 
ure [5] compares experimental and calculated dephasing 
rates which show comparable monotonic increases with 
increasing temperature over this range. 



IV. CONCLUSION 

Approximate quantum dynamics methods which are 
based on combining path integral expressions for forward 
and backward propagators and linearizing in the differ- 
ence between these paths provide a popular approach for 
incorporating some quantum effects with in a trajectory 
based framework. In this paper this type of approach has 
been implemented to study evolution of the off-diagonal 
elements of the vibrational density matrix which describe 
vibrational pure dephasing. In this approach trajecto- 
ries of "bath" variables that influence the quantum sub- 
system evolve on the "mean potential" associated with 
the relevant states. This quasi-classical mean potential 
dynamics emerges naturally from the linearization the- 
ory which also gives a unique prescription for the initial 
conditions of these trajectories. With this approach the 
Wigner transform of the full quantum equilibrium den- 
sity operator is found to provide the initial phase space 
distribution for the quasi-classical trajectories. In this 
paper we have explored the accuracy of various methods 
for obtaining approximations to this quantum initial den- 
sity. Our measure of quality here is the reliability with 
which these different methods can reproduce vibrational 
pure dephasing dynamics over a range of temperatures. 
Thus, for example, we found that the classical initial 
phase space distribution gives pure dephasing rates that 
agree with experiments in higher temperature crystals. 
As a consistency check the FK-Wigner approach for ap- 
proximately sampling the full quantum phase space dis- 
tribution which is developed in detail in this paper gives 
results that also reproduce the experimental findings in 
these higher temperature solids. However, in low temper- 
ature solids we find that the classical initial phase space 
distribution combined with the quasi-classical dynamics 
can not reproduce the experimental dephasing dynam- 
ics. When the quasi-classical trajectory initial condi- 
tions, however, are sampled from the FK-Wigner approx- 
imate quantum equilibrium density excellent agreement 
between experimental and calculated dephasing rates is 
found. 

There has been considerable theoretical work exploring 
vibrational dephasing in condensed phases. For example, 
Zewail and co-workers give a detailed exposition of the 
perturbation theory results applied to the anharmonic 
oscillator in a harmonic bath modeli^i^. Skinner and co- 
workers explored dephasing for a model two level system 
coupled to a harmonic bath where they obtained non- 
perturbative results^^. Generally when the Debye-like 
spectral density is employed in these different theoreti- 
cal results the dephasing rates they predict have a very 
strong temperature dependence with 7 ~ T 7 . Similar 
sensitivity of dephasing rates to spectral density has been 
found in numerical simulations employing the stochastic 
classical trajectory approach 69 ' which ignores quantum 
effects. As discussed in section Mil our realistic micro- 
scopic model calculations, which make no assumptions 
about interaction strength, anharmonicity, or assump- 
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FIG. 5: Comparison of experimental and calculated vibrational pure dephasing rates as functions of vibrational quantum state. 
Each panel compares FK-Wigner and Classical initial condition sampling calculations with experimental results. The different 
panels report results for various temperatures: left panel, T—2.6K; center panel, T=20K; right panel, T=32K 





FIG. 6: Left panel compares pure solvent (solid krypton) pair distribution functions at T=32K computed using Monte Carlo 
sampling with the Classical Boltzmann distribution (solid curve), and the FK approximation to the Wigner transformed 
equilibrium distribution (dashed curve). Right panel presents the same results but for T=2.6K. 
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FK-Wigner T=2.6k 
FK-Wigner T=20k 
FK-Wigner T=32k 
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FIG. 7: Pair distribution functions for solid krypton at various 
temperatures computed use FK-Wigner approximate quan- 
tum equilibrium distribution. 



tions about the nature of the underlying spectral density, 
and include quantum nuclear initial conditions effects, re- 
produce the very weak temperature dependence observed 
in the experiments. Increasing the temperature by about 
a factor of 10 (2.6K to 32K) results in a modest increase 
in dephasing by about a factor of two for most vibra- 
tional levels. Using their non-perturbative theory and a 
pseudo-local mode spectral density model which assumes 
that the vibrational dephasing occurs by coupling to low 
frequency modes that arise from, for example, hindered 
rotation and translation of the vibrating impurity in the 
crystal, Skinner and co-workers (and other o 71 i 72 ' 73 i 74 ) ob- 
tained various Arrhenius like forms for the temperature 
dependence of the dephasing rate. Such forms are con- 
sistent with the much weaker temperature dependence 
that we find in our calculations and in the experiments. 
A "back of the envelope" estimate using a single Arrhe- 
nius form, and assuming a temperature independent pre- 
exponential factor gives an activation energy or pseudo- 
local mode frequency of about 20 cm -1 . In future work 
we will explore the nature of these pseudo-local modes re- 
sponsible for this dephasing dynamics by examining the 
microscopic motions responsible for dephasing in our re- 
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FIG. 8: Comparison of experimental (left panel) and calculated (right panel) vibrational dephasing rates as functions of initial 
vibrational state superposition (0, v) for iodine in solid krypton. 



alistic simulation model^. 

The study of dephasing dynamics in condensed phase 
systems is just one example that highlights the impor- 
tance of reliable methods for including quantum effects 
in condensed phase calculations. With in the frame- 
work of linearization these quantum effects are incorpo- 
rated principally in the initial distribution sampling. At 
present there are no direct ways to generate the Wigner 
transform of the thermal density operator and sample 
the exact full quantum distribution. Thus all current 
methods for implementing linearized dynamics for gen- 
eral applications use local Harmonic approximations. As 
shown in section Iff Bl the FK- Wigner approach employs 
a thermodynamic variational calculation to parameterize 
a locally quadratic approximate form for the Euclidean 
action appearing in the partition function. This form is 
then employed to approximate the Wigner transform of 
the density operator. With such an approach no local 
harmonic approximation is made to the potential. Thus 
the method incorporates global information from the par- 
tition function integral in order to fit the gaussian ap- 
proximate form for the density. This leads to a form for 
the density operator that contains frequencies computed 
from a gaussian smearing of the interactions which can be 
implemented in a very computationally efficient means. 

This situation is contrasted with that of alternative 
methods such as that of Shi and Geva^i, in which a 
local Harmonic approximation to the bare potential is 
made. While these methods can use full PIMC calcu- 
lations to accurately represent diagonal elements of the 
density matrix, using a local harmonic approximation for 



the off-diagonal elements can degrade the results. Thus 
the studies outlined in this paper suggest that the FK- 
Wigner approach is the method of choice for a balance 
of computational efficiency and accuracy. 
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APPENDIX A: COMPUTATION OF POTENTIAL 
AVERAGES OVER TRIAL VARIATIONAL 
DISTRIBUTION 

In order to apply the variational result in Eq. (|2"0|) the 
various potentials averaged over the trial density defined 
in Eqs. (f2"Tj) - (|2"5|) must be computed. These calculations 
are most conveniently performed in Fourier space, thus 
we write 

V(R(0)) = V(R) = J j^y d V{ky k R (Al) 
Further since R(t) = M _1 / 2 U?7(t) we find that 



drj c -rl( Vc ) TfTJ I dRerj nj f dlmr] n 



xe -/3P^+^ 2 (')c)][(Rc t; „ 3 ) 2 + (Im^„ J ) 2 ] f°° dk y( fc yfe r fl Ce 2E J -6 J -(Rs%»j) ( A2 ) 
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with bj = i(k T Wl x / 2 XJ)j. The integrals over the real We can compute the Gaussian smeared local harmonic 



and imaginary parts of the r] n j can be performed analyt- 
ically using the same result that lead to the expression 
in Eq. l|2"Tj) , and the final result has the form 



approximate trial potential in a similar fashion obtaining 



Vi Hal (R c ) = / 

J — i 



dR 



-e 2 



\(R-R e ) T AT 1 {R )(R-R a ) 



(V(R)) 



1 



trial 



Zfrial 
d 



dR c 



X 



{2-Kh 2 f3) d l 2 

phwj{R^/2 
smh.phjJj{R c )/2 

dk 



|27rA|V2 

l(R - R C ) T M 1 / 2 K(R C )M 1 ^(R - R c ) + L(R C 



n 



lj2 M l /2K ^M- /2 A ij+ L(R c ) 



(A10) 



(2tt)' 



jV(k)e l 



= e 2 



(A3) 



where we have used the following 



APPENDIX B: VARIATIONAL CALCULATION 
DETAILS 

To simplify notation we define 

d 



Pp(R c 



l 



fjhw 3 {R c )j2 
(2nh 2 /3) d / 2 sinh/3^(i? c )/2 



n 



i 



fc T Afc 



defining the smearing width matrix A as 
A = M- 1/2 UAU T ]Vr 1/2 

with 

2 



(A) y =E 



(A4) 



(A5) 



(A6) 



and 



X&c) = lY, M } / * K * A U M : 



1/2 



so we can write 



Z tna i = J dR c e-^ R ^P p {R c 



(Bl) 



(B2) 



(B3) 



Since V(k) = FT{V} and e -§fc T A(fl c )fc 
FT{|27rA|- 1 / 2 exp[-iJ2 T A- :L ^]}, and defining Va(# c ) 
as follows, we see that 



V A (R C ) = 



dk 



-co (2tt 



ik T Ro e ~^k T A(R c )k _ 



f 1 1 

FT' 1 i FT{V} * FT{\2ttA\~ 1/2 exp[- -R t A' 1 R}} \ (R c ) 



Following Feynman and Kleinert^ we proceed to opti- 
mize the right hand side of Eq. (f2"0|) , / = exp[— i(S — 

Striai)triai]Z t riai = e~ //Zt " a! Z tr i a i , with respect to the 
variational parameter functions K(_R C ) and L(R C ). Here 
we have defined 

I = (3 JdR c [V A (R c ) - x(Rc) - L(R c )}e-P L{R ^P p (R c ) 

(B4) 

Partial functional differentiation of /[K, L\ first with re- 
spect to L (fixed K) yields 

5 L f = e- T ' z ^' [{I/Ztrial + l)S L Z t riai - 5 L I] (B5) 



(A7) 



and since, for example, 



Using the Fourier convolution theorem (FT {fx * /^j = 
FT{fx}FT{f 2 }, where (/i*/ 2 )(<?) = / dy fx(y) f 2 (q - y)) 
we can write 



5,,I 



dR r 



81 



SL(R C 



-(R C )SL(R C ) 



(B6) 



(V(R)) 



1 



trial 



Ztrial 
d 



dR c 



(2nh 2 f3) d / 2 



-/3L(R C 



n Phw 3 (R c )/2 
X l) iS inhf3h, 3 (R c )/2 VA{Rc) (A8) 

so Va(R c ), defined above, can be interpreted as the Gaus- 
sian smeared potential since 



V A (Rc) = 



dR 



\2-kA\ 1 / 2 



3 -i(iJ-fl. c ) T A- 1 (_R c )(_R-_R c 



V(R) 
(A9) 



Where the partial derivative notation means regard the 
functional I[L(R C )] = I(L(Rl), L(R 2 C ), L(R™)) in the 
limit n -> oo as a multidimensional function of L evalu- 
ated at various points R\, and differentiate with respect 
to L(R C ) at the point R c while keeping all other values 
of L(R' C ^ R c ), and the values of the functions K at all 
points fixed. Thus we find 

S L f = e-'l^if J dR c [V A (R c )- X (Rc)-L(R c )} 
e-^ R ^P (R c )5L(R c ) 

-131 /Z tri ai J dR c e-P L ^P p (R c )5L(R c )} (B7) 
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The only way to get both terms in the curly brackets 
in this expression for S^f to vanish for arbitrary 8L{R C ) 
is if the quantity [Va(R c ) — x{Rc) — L(R C )] is zero for 
each R Cl giving the following local relationship between 
the parameter functions characterizing the extrema with 
respect to L 



L e (Rc) = V A (R C ) -^J2 M i /2K ij( R -) A ^( R -) M : 



,(B8) 

With this result 1 = 0, and we now must optimize 
/[K, L e ] = Z tr iai[K., L e ] with respect to variations in the 
functions K with L fixed at L e . Now 



J trial 



!K1 



dR c 



(B9) 



and Wtriai for L = L e is obtained by substituting Eq. (|B8|) 
into Eq. (|28|) to give the following expression for the trial 
effective potential that depends only on the local curva- 
ture parameter matrix K 



Moreover, since 



1 2x^ 1 
cotn7rx = 1 > —~ — 



n—1 



(B15) 



it is easily shown that 



l_ { ^ coth ^_ i} 



(B16) 

and thus the first term on the right hand side of 
Eq. (|B12[ ) vanishes identically using the results around 
Eq. (|A5jl . Thus Eq. (|B12p gives the optimization condi- 
tion as (dWtriai/dAij)n = 0, since in general dAij/dKij 
is nonzero. 

Differentiating Wtriai with respect to a particular ele- 
ment of A keeping all others fixed (and pretending that K 
can be held fixed during this process consistent with the 
application of the chain rule) we find, using the Fourier 
space form of Va in Eq. (|A7[) . that 



W tr i a i(R c ) 



f3hw 3 {R c )/2 



(BIO) 



d 

-- Vln, 

Pj^i \smh/3haj 3 (R c )/2 

+v A (Rc) -\Y. m I /2k ^Rc)a 13 \r c )m. 



trial 



1 „ r dk 



0A fk , 2 L J (2tt) 



-M)! 2 K yk ,Ml! 2 } 



k r A(fl c )k 



(B17) 



Functional optimization with respect to each of the func- 
tions Kij requires 

5 Ki ,Z trial = J dR c ^^- ) SK l3 (R c ) 

SWt, 



(Bll) 







For this to be true for arbitrary variation functions 
SKij(R c ), (dW t riai/dKij)(R c ) must vanish locally for 
each function Kij. Since A is a complicated function of 
K given by Eq.jM} we write Wtriai = W / t ria ;(K, A(K)) 
and apply the chain rule to obtain 



dW t 



trial 



dKi 



dWtrv 

dKa 



dWtr 



dA { 



(B12) 



After some straightforward manipulations we find that 

' dWtrial ' 



dK 3 , k , 



l -MfAj, k ,Ml!> 



y^ { ^L coth 1 (bis) 

1 9 9 Mir.- . I 



1=1 



dK 



j'k' J A 



Further since lo 2 = U T KU is diagonal, so u>i = 
Ej E fc U^KjkUu] 1 / 2 and thus 



9Kj'k> J A 2w, 



U^U m 



(B14) 



Packing all these terms into a dxd matrix, the extremum 
condition requires that all elements are zero and the re- 
sult can be written in matrix form as follows 



rfk 



(27T) d 

It is easy to show that 



* kT ^(k)kk T e^ kTA ^ k = M 1 /2 K M 1 / 2 

(B18) 



FT 



<hV X = V(k){ik){ik) T (B19) 



dRdR T J 



So using the same manipulations based on the convolu- 
tion theorem that lead to the result in Eq. (|A9[) we obtain 



K(R C )=J ^-^D(i?)e-^-^) TA -V.)(^) 

(B20) 

Where D(i?) is the mass weighted Hessian 

D(i?) = M- 1 / 2 (g^^j (R)M-^ (B21) 



APPENDIX C: GENERAL APPROACH FOR 
GAUSSIAN SMEARING PAIR POTENTIALS 



The left hand side of Eq. (|B18p is conveniently evalu- 
ated in Fourier space if the interaction potential is pair- 
wise additive V(R) = v%3 ( r u )' an( i eacn P a i r func- 
tion of the inter-particle distance r*y = \Ri — R 3 \ is fit to 
a sum of Gaussians v l i{n 3 ) = Y^ k a k ex P( — \b k 3 r 2 ). In 
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k-space this potential has the form 



3/2 



Using this result we can rewrite the general (mi, ny) ma- 
trix element of Eq. (|BT8]) giving 



xe 



-\k,\ 2 /2bl 



(CI) 



(flcm-Hc„)l. I, 



-ik,i 2 /2fcr 



3 3 



i=l j=l 



r 



Constructing the symmetric 3x3 matrix 7"™ for particles 
n and m with elements of the form 



A, 



[A 



mn 

mxnx ~r ^nsmsj ~r J-/ 



(C2) 



and 

Taiy — \A-mx7ny H - A nX ny] \A-mxny A nX my] (C3) 



the remaining k m integral can be readily performed yield- 
ing the following expression for the 3x3 K mn sub-matrix 



K mn = {M m M n )- l ' 2 Y - ° fc = 
^y/(br) 3 detj 

X exp[-i(i? cm - i? c „) T (7" m ) _1 (^cm ~ Ren)} 



[(j mn )~ 1 — u m ™( u m ™) T ] 



r 



(C4) 



Where u" m = (7 mn )~ 1 (-R cm - Rcn). This approach is 
readily simplified to compute Va(R c )- Proceeding in a 
similar fashion to that outlined above employing the k- 
space expression of Eq. l|A7|) we obtain 



APPENDIX D: LOCAL HARMONIC 
APPROXIMATE DENSITY OPERATOR AND 
ITS WIGNER TRANSFORMATION 



v A {R c ) = yy^^= 

X exp[-l(i? cm - RcnfiT^-^Rcm - Rcn)] 



(C5) 



Thus we define the density operator with the varia- 
tionally optimized trial harmonic form as 



-pH 
trial 



q{Ph)=q" i r ph 

dq' I dq"\q')(q"\ I dq c I V[q(T)]5(± / q(r)dr - q c ) 

J Jq(0)=q' P n JO 

0H 



r 1 

x exp[— — 



1 



1 



dr{-q(ry q(r) + L(q c ) + -(q(r) - q c Y K(q c )(q(r) - q c )}} 



(Dl) 



or writing the (^-function in integral representation we 
have 
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0H I dq' dq"\q')(q"\ dq c ^ V[q{r)] (D2) 



''trial 



<>M- T I dr{- q (rf q (r) + L(q c ) + -(q(r) q c ) T K(q c )(q(r) - q c ) + (pk T (q(r) q c )}] 



I 

Transforming to normal mode vectors rj = XJ T q where U matrix so that K(q c ) = Uw 2 (<7 c )U T as in Eq. p5|) , and 
diagonalizes the variationally optimized local curvature defining y(r) = (t](t) — r) c ) this expression becomes 



e 



= Idrf I drf'\rf)(rf'\ I d Vc e-^) I / v[y[T)] 



I rid 



(2tt)° 



V(0) = (ri' -Vc) 



1 f l3h 1 2i 
x exp[-- J dT-{y(r) T y(r) + y{T) T uj 2 y( T ) + (-) K T y(r)}] (D3) 



where k = XJ T k. However, the path integral 

rx{(3h)—x pf3h ^ 



1=1 V[x(r)]exp[-^ I dT^-{i 2 (T) +uj 2 x 2 {t) +2cix(t)}} 

x(0)=x Jo 2 



can be readily evaluated yeilding 



_ / OJ \ 1/2 , lo(x' 2 +x 2 ) ujxx' 

V27rfi^mhpfiw) ) GXP ^~ 2ntanh(/3fiw) + ftsinh(/3ftw) J ^ ' 

2 

x exp[-^{(/3W2) -tanh(/3fiw/2)} + ^-(a; + x')tanh(/3^/2)] 



I 

With a = 2iK n /P as in Eq. (|D3p , I has a Gaussian form in Eq. (|D3p can be performed analytically. The final result 
K n and the Gaussian integrals over these components in of these manipulations is 



Here we define the quantity 

a, = coth /3ftwi/2 - 2/phui (D6) 

The above locally Gaussian form in the normal mode 
variables rj' , and 77" for the Boltzmann operator is partic- 
ularly convenient for approximately evaluating thermally 
averaged time correlation functions in the linearized 
approximatio n 14 ' 20 ' 41 where quantities like Wigner trans- 
forms of products of operators of interest with the Boltz- 



mann operator i.e. 

(P0O)w(Q,P) = J dz{Q + z/2\ppd\Q-z/2)e~* Pz 

(D7) 

are required. For simplicity we will focus on the Wigner 
transform of the Boltzmann operator itself. In the 
applications to computing linearized approximations to 
time correlation functions mentioned above the quan- 
tity {p[}0)w{Qt P) serves both to provide a distribution 
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for initial conditions for classical-like trajectories and in- 
cludes the measurement of the initial operator. Depend- 
ing on the nature of the initial operator O, the prod- 
uct Wigner transform may be dominated by the Wigner 
transform of the Boltzmann operator. This is the case 
in the applications to vibrational dephasing that we will 
report in a subsequent publication 7 - as the operator in 
question just involves the vibrational states of the so- 
lute chromophore. The solvent initial degrees of freedom 



must thus be sampled from the Wigner transform of the 
solvent Boltzmann operator. Here we thus explore the 
accuracy of the Feynman-Kleinert approximation for the 
initial Wigner phase space distribution and the underly- 
ing approximation to the thermal density matrix. The 
Wigner transform of the trial Boltzmann operator writ- 
ten in the local normal mode phase space representation 
thus has the form 



= j (D8) 

xexp[ a t n {Ql Vci > Jcxp[ nut Qlli 

I 
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